g L / ^ vy 5 


Theoret. Comput. Fluid Dynamics (200D M: 399-422 Theoretical and Computational 

Fluid Dynamics 

© Springer- Vcrlug 2(X)1 


Turbulence Modeling Effects on the Prediction of 
Equilibrium States of Buoyant Shear Flows* 

C.Y. Zhao and R.M.C. So 

Department of Mechanical Engineering, The Hong Kong Polytechnic University, 

Hung Horn, Kowloon, Hong Kong, P.R. China 


T.B. Gatski 

Computational Modeling and Simulation Branch, 
NASA Langley Research Center, 
Hampton, VA 23681-2199, U.S.A. 


Communicated by M.Y. Hussaini 


Received 21 April 2000 and accepted 21 February 2001 


Abstract. The effects of turbulence modeling on the prediction of equilibrium states of turbulent buoyant 
shear flows were investigated. The velocity field models used include a two-equation closure, a Reynoids- 
stress closure assuming two different pressure-strain models and three different dissipation rate tensor 
models. As for the thermal field closure models, two different pressure-scrambling models and nine dif- 
ferent temperature variance dissipation rate (g«) equations were considered. The emphasis of this paper 
is focused on the effects of the ^-equation, of the dissipation rate models, of the pressure-strain models 
and of the pressure-scrambling models on the prediction of the approach to equilibrium turbulence. Equi- 
librium turbulence is defined by the time rate of change of the scaled Reynolds stress anisotropic tensor 
and heat flux vector becoming zero. These conditions lead to the equilibrium state parameters, given by 
P/ E> p#/eu, R - (0*/2£o)/(k/s), Sk/e and G/e , becoming constant. Here, P and Pn are the production 
of turbulent kinetic energy k and temperature variance 0 2 , respectively, e and e h are their respective dis- 
sipation rates, R is the mixed time scale ratio, G is the buoyant production of k and S is the mean shear 
gradient. Calculations show that the ^-equation has a significant effect on the prediction of the approach 
to equilibrium turbulence. For a particular ^-equation, all velocity closure models considered give an 
equilibrium state if anisotropic dissipation is accounted for in one form or another in the dissipation rate 
tensor or in the /--equation. It is further found that the models considered for the pressure-strain tensor 
and the pressure-scrambling vector have little or no effect on the prediction of the approach to equilibrium 
turbulence. 


Introduction 

Homogeneous turbulent flows played a central role in the modeling and analysis of complex inhomogeneous 
turbulent flows (Rogallo and Moin, 1984; Speziale. 1991). The reason is that these homogeneous flows pro- 
vide a great deal of insight into key parameters characterizing turbulence in a simplified setting which quite 
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often gives rise to closed-form solutions in the commonly used turbulence models. For turbulent shear flows 
without buoyancy, Lee and Reynolds (1985) and Rogers et ai (1986) have conducted direct numerical sim- 
ulations (DNS) of such flows, and experimental investigations were carried out by Harris et ai (1977), Rohr 
et aL ( 1988a) and Tavoularis and Karnik ( 1 989). These investigations provided much-needed insight into the 
understanding of the evolution of homogeneous turbulence and its approach to equilibrium. 

For turbulent buoyant shear flows, DNS have been performed by Gerz et ai ( 1 989), Gerz and Schumann 
(1991) and Kaltenbach et ai ( 1 994), while Piccirillo and van Atta ( 1 997) and Rohr et ai ( 1 988b) have exper- 
imentally studied the buoyancy effects on stratified turbulent shear flows. According to these investigations, 
stable stratification weakens isotropization. Therefore, the isotropic dissipation rate model is not adequate 
to describe the evolution of buoyant homogeneous turbulence. Recently, So et ai (1999) conducted a nu- 
merical simulation using different turbulence models and assessed the performance of certain anisotropic 
dissipation rate models. Their investigations revealed that there were difficulties still in the prediction of 
counter-gradient heat flux and the onset of internal gravity waves in buoyant shear flows. However, it should 
be pointed out that all the investigations mentioned above were only carried out for the short time period (or 
the near-field region of experiments), therefore, turbulence evolution for the long time was not reported due 
to experimental and numerical difficulties. 

An important property of homogeneous turbulent shear flows is the appearance of the dynamic state pa- 
rameters, which tend to approach equilibrium values in the long time limit. The equilibrium states provide 
an important benchmark in the calibration of closure models. For non-buoyant turbulent shear flows, Abid 
and Speziale ( 1993) calculated the equilibrium states for channel flows and homogeneous turbulent shear 
flows using Reynold s-stress closures. The fixed points associated with the equilibrium states for several ho- 
mogeneous shear flows were determined from bifurcation diagrams (Speziale and Mac Giolla Mhuiris, 1 989; 
Speziale, 1991; Speziale eta/., 1996), and they were used to assess the stability of higher-order models 
and their ability to predict the correct equilibrium values. Recently, using representation theory, Jongen and 
Gatski (1998) also showed that equilibrium turbulence of homogeneous shear flows is defined by the state 
parameters. P/e and Sk/e. Here, P is the production of turbulent kinetic energy k and S is the mean shear 
rate of the mean flow 7 . They further derived a general algebraic relation between the state parameters based 
solely on the form of the pressure-strain rate model, and without having to specify a modeled ^-equation 
lor the dissipation rate of k. F heir analysis gave the same equilibrium values as those predicted by the 
Reynolds-stress turbulence closures. 

Tavoularis and Corrsin (1981) conjectured that equilibrium states also exist for buoyant shear flows, how ; - 
ever, they take a longer time (or distance) to achieve. Due to the interaction between shear and buoyancy, the 
approach to equilibrium turbulence tor buoyant shear flows is much more complicated. Up to now, even the 
parameters that characterize the equilibrium states are not known precisely and there was a lack of reliable 
experimental data on the equilibrium states of turbulent buoyant shear flows. 

Numerical modeling of buoyant shear flows, even under the assumption of incompressibility, is quite 
a bit more complicated than pure shear flows. In addition to closure models for the Reynolds-stress 
equation, other models for the equations governing the transport of the Reynolds heat-flux and the asso- 
ciated dissipation rate (£#) of the temperature variance (O 2 ) need to be invoked (Launder, 1978, 1989). 
Furthermore, buoyancy effects should also be modeled into the pressure-strain term based on the argu- 
ment that buoyant forces enter into the Poisson equation for the fluctuating pressure. A modeled buoy- 
ancy term was also added to the ^-equation to account for buoyancy effects in the calculation of e. 

It should be noted that this treatment of buoyancy effects on the turbulence field is rather incomplete, 
because the interaction between shear and buoyancy is highly nonlinear. This treatment also ignores 
the tight coupling between the Reynolds-stress equation and the heat-flux equation, because it simpli- 
fies the coupling between the different terms in these equations. Since the exact form of the modeled 
^-equation is not important in the prediction of the equilibrium states of pure shear flows, would the 
same conclusion be true lor the -equation for buoyant shear flows. In view of this, there is a need 
to explore the ability of commonly used turbulence models to predict the equilibrium states of buoyant 
shear flows. 

An analysis along the lines proposed by Jongen and Gatski ( 1998) was used to investigate the equilibrium 
state parameters of buoyant turbulent shear flows (Zhao et ai, 2001). The requirements that dbjj/dt = 0 
and dr/, yd t — 0 l ead to P/e , Sk/e , Po/e.o, R = (6 2 /2eu)/{k/e) and G/e becoming constant at equilibrium. 
Here, bq = itjUj/lk — 8,]/ 3 is the scaled Reynolds stress anisotropic tensor, q, = u~jb/(k0 2 ) ]/2 is the scaled 
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heat flux vector, m,-m/ is the kinematic, Reynolds stress tensor, u,0 is the heat flux vector, Ph is the pro 
duction of the temperature variance 0 2 , % is its dissipation rate, R is the mixed time scale ratio and G 
is the buoyant production of k. Therefore, the objective of this paper is to examine the effects ot turbu- 
lence modeling on the prediction of these equilibrium states. It can be classified into the following tasks. 
The first task is to investigate the effect of the modeled /^-equation on the prediction of equilibrium states. 
A second task is to study the effect of modeling the pressure-strain tensor. The third task is to examine 
the effect of the model invoked for the pressure-scrambling vector. Finally, a tourth task is to investigate 
the influence of anisotropic dissipation modeling. Since this model affects the modeled f -equation, the 
effect of this equation in the context of anisotropic dissipation modeling will be examined also. Two dif- 
ferent turbulence closure schemes are used in this study, a k-e and a Reynolds-stress closure. The k-e 
closure is used to study in detail the effect of the modeled ^-equation, while the Reynolds-stress closure 
is used to examine the effects of pressure-strain models, pressure-scrambling models and dissipation rate 
tensor models. 


The Modeled Equations 
Reynolds Averaged Equations 

The flow considered is a homogeneous buoyant turbulent shear flow' where the Boussinesq approximation 
is assumed valid. According to the studies of Sommer and co-workers (1995, 1997), two-equation heat flux 
models are not appropriate for buoyant shear flows, because they failed to predict the onset of counter- 
gradient heat flux and the presence of internal gravity waves correctly. Therefore, a Reynolds heat-flux model 
has to be assumed instead. Invoking the Boussinesq approximation, the incompressible Reynolds-stress and 
heat-flux modeled transport equations can be written for a homogeneous shear flow as 


d iijitj 

a r~ 

d iiiO 

~dT 


p )j -r n ij — s ij + g ij , 

(K~) dUi 

— U i U j UjO — h + O Oi « 

J (kj ax j 


( 1 ) 


( 2 ) 


where U\ and «,■ are the / th component of the mean and fluctuating velocity, <r) and 0 are t he mean and 
fluctuating temperature, t is time, ,v, is the /th component of the coordinate, P// = -«/«* (3Uj/dx k ) - 
Tfiuk[dUi/dxk) is the shear stress production tensor, C,;/ = — figiUjO - figjiiiO is the buoyant production ot 
the Reynolds stresses and Gm = is the buoyant production of the heat fluxes. Here, is the coeffi- 

cient of thermal expansion of the fluid and Hi is the gravitational vector. The pressure-strain tensor n y and 
the dissipation rate tensor in (1 ) need modeling. The <t>w term includes both the pressure-scrambling vec- 
tor and the molecular dissipation of the heat fluxes. Its writing was suggested by Launder (1978). Therefore, 
its modeling should include the effect of molecular dissipation on the heat fluxes. These models are briefly 
discussed below. 


Pressure-Strain Models 

(S ) 

For buoyant turbulent flows, the term n,y can be represented by three parts; they are the slow' part n,. . 

the rapid part n) R) and the buoyant part n^’. The slow and rapid parts can be modeled as in non-buoyant 
flows and the buoyant part can be modeled by relating it to the generation/destruction of the Reynolds 
stresses by buoyancy (Launder, 1978). Two different models for nJ ; S) + n' ; R) are investigated; they are the 
linear IP model of Launder etal. (1975) and the nonlinear SSG model of Spezialeem/. ( 1991 )• These models 
plus the buoyant model, n|?’ = C 6 (C,; - 2C5„/3), can be written as follows. The complete expression for 
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invoking the IP model for nj?’ + n| R) , is given by 

n v = _ C\ebjj + CikSjj + Cyk ^b, k S jk +b jk S, k - -b m „S mn 8jj^ 

+ C 4 k {b ik W jk + b jk W ik ) - C 6 (cij - . (3) 

where C\ = 3.6. C2 = 0.8, C\ = 1 .2, C4 = 1 .2, C(, = 0.3 and G = G„/2 is the buoyant production of k. The 
corresponding expression for n y . invoking the SSG model for njj 1 + n : R) , is 

n v = - (o + C* P) ebjj + C 2 kSij + C 4 k (b ik S Jk + b jk S lk - ~b mil S,„„8,j"j 

-C A k (b ik W jk - b jk W ik ) + C 5 s (b ik b kj - l -b mn b nm 8ij^ - C 6 ^Gy - ^ GS^ . (4) 

The model constants are specified as C\ = 3.4, C 2 = 0.36, C 3 = 1 .25, C 4 = 0.40. C 5 = 4.2, C* = 1 .8 and 
C 6 = 0.3. in these equations the kinematic strain rate (S/y) and rotation rate ( W l} ) tensors are defined as 
Sjj = (dUi/dxj + Wj/dxi) /2 and Wjj = [dUjfdXj — dUj/dXj) /2, respectively. 

Dissipation Rate Models 

The dissipation rate tensor £,y can be modeled using three different approaches. They are the isotropic model 
(Kolmogorov, 1941) and two different algebraic models (Hallback et al. , 1990; Speziale and Gatski, 1997). 
These models are given here for reference without derivation. The Kolmogorov (1941) isotropic model is 
given by 


2 

3 eS U ■ 


(5) 


while the anisotropic algebraic models of Hallback et al (1990) and Speziale and Gatski (1997), designated 
as HGJ and SG, respectively, can be written in the general form as 


£// — “£<$// -\-2ed l 


V * 


(6) 


In the HGJ model, djj is expressed as a series of /?,y , while in the SG model, d,j is expressed in terms of Sa 
and Wjj . For the HGJ model, the expression for djj is given by 


da = - 


* G n _ 5) hij ~ a ' ( h,khk J ~ ^ - 


where a' = j has been assumed. The SG modeled equation for djj is given by 


dij ~ -d’fif 


k 

-s„ 


C FS + P/e- I | 


-2(Si k W k j+S jk W ki ) 


+2 C fie 


^2 (^Sj k S k j - -S k iS k i8ij^j . 


C f 5 + P/e - 1 ] e* 

where a 3 - 0.6. A\ = -pj-ar.3 + and A 2 — jjff3 - f j-. The coefficient C, u is defined as 


(7) 


( 8 ) 


C klt = 


15 (C f -5 + P/s — I ) 


1+2 


C f 5 + P/e — I 


) ^ _ 5(c,5+^/5-l) 


(9) 


Two new constants are introduced and they are C f5 = 5.80, A 3 = {^<*3 - 2,. The dimensionless strain invari 
ants are defined by m = {S i jS ij ) l/2 (k/e) and G = (W ij W ij )^ 2 (k/e). 
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The Sjj model requires knowledge of s for closure. For the HGJ model, the e-equation solved is the 
conventional equation, which after modification for buoyant flows becomes 


de 

dt 


2 

C«| 7 P - c, 2 - Crf jPgwB • 

k k k 


( 10 ) 


where P is one-half the trace of Pjj. The model constants take on the following values: C?\ = 1.44, 
C e0 = 1 .83 and C& — 1 .44. Note that (10) does not include the effect of djj , whose influence only contributes 
to the calculation of Sjj through (6). This, however, is not the case for the SG model, and the E-equation 
solved is modified to give 

2 

— = C h \ 7 P - 2( 1 + y)edjjSij - C,i E — - , ( 1 1 ) 

dt k k k 


where y = | (n«3 - jf) and the constants take on values C*i = 1.0, C f 2 = l- 8 ^ and 0-3 = 1-44. ln th ' s 
equation the effects of djj are accounted for through the term involving both dy and S, r Thus, the anisotropic 
effects of £ij are given by (6) as well as accounted for in the determination of e. This is a major difference 
between the HGJ and SG models. 


Heat Flux Models 

Just like riy, the term can again be partitioned into three parts, a slow part 4>^\ a rapid part O,#, and 
a buoyant part Usually, the Monin (1965) model is invoked for However. Lumley and Khajeh- 
Nouri (1974) proposed a quadratic modification to Monin’s (1965) model to account tor the anisotropy of 
the heat flux. Their proposal leads to = —C\n£Uj9/k ■ 2C'- lf y:h ll u l 0 / k. The rapid part is usually mod- 
eled by = b'^dUk/dXm) with different models proposed for Recently, Craft and Launder ( 1989) 
derived a fairly complete expression for b'^ based on the realizability condition in teims of Schwarz s in- 
equality for the heat fluxes and by imposing additional conditions. Their expression involves quadratic terms 
of bjj in addition to terms linear in the heat flux. If all quadratic terms involving turbulent stresses and heat 
fluxes are neglected and the model 4*^ ' = (Launder, 1978) is invoked for the buoyant part, the 

complete model for <t>« can be written as (So and Speziale, 1999) 

<t> w = —C\nji(iO + C2nitfi-r— + CyuitjO—d- + CsvgjfiO 2 , ( 12 ) 

k oXj oXi 

where C 50 = Ciu is recommended by Launder (1978). If the term involving C 30 is set to zero as suggested 
by Launder (1978), the other constants would take on the following values, C w = 3.28 and C 29 = 0.40. 
On the other hand, including the Cw term leads to the recommended constants C\o- 3.28, C 20 = 0.80 and 
C if) = -0.20. 

The linear model given in (12) is by no means unique. For example, Shabany and Durbin (1997) sug- 
gested a slightly different linear model for O ^ 1 . Their proposal consists of three terms, the Cy t and Cvt term 
plus a term involving the mean temperature gradient, CwUjUj(d&/dxj). However, in their calculations, they 
have set C<w - 0. Therefore, in reality, their suggested linear model for is the same as that given in ( 1 2). 
A comparison of the effects of the different terms in (12) including a C49 term on heat flux modeling has 
been attempted by Wikstrom et al. (2000). Their comparisons were made with DNS data, whenever avail- 
able, obtained for a homogeneous shear flow, a plane channel flow and a heated cylinder wake. In addition, 
they have also considered models where the Cj« term was modified by the time scale ratio R. Their results 
showed that if the C 49 term is included in (12), then C 2 n and C 3 # have to be set equal to zero. Otherwise, 
the model performance would deteriorate significantly. A model that performs well for all test cases was one 
given by modifying C\» with R and setting all other constants to zero. Other models perform differently for 
the test cases considered. 

As a first attempt, the present study focuses on ( 1 2) only. It should be pointed out that the proposed model 
includes the effect of molecular dissipation on the heat fluxes and is consistent with those used by other re- 
searchers. All calculations reported in this paper, unless specified explicitly, were carried out assuming C.w 
to be zero. The effect of the C 20 term on the prediction of equilibrium states is investigated separately. 
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Table 1 . A summary of the ^-equation model constants proposed by various researchers. 


M ode led -eq uat ion 

0/1 

0/2 

0/3 

0/4 

O/s 

NLL (Newman et ai, 1981) 

2.00 

0 

0 

2.02 

0,88 

LL (Elghobashi and Launder. 1983) 

1.80 

0 

0 

2.20 

0.80 

GJK (Gibson et ai, 1987) 

0 

1.70 

1.40 

2.00 

1.80 

z\Y (Yoshizawa, 1988) 

1.20 

0 

0.52 

1.20 

0.52 

JM (Jones and Musonge. 1988) 

0 

1.70 

1.40 

2.00 

0.52 

NK (Nagano and Kim, 1988) 

1.80 

0 

0.72 

2.20 

0.80 

C’L (Craft and Launder, 1989) 

0 

1.30 

0 

2.00 

0.92 

YNT ( Youssef et ai, 1992) 

1.70 

0 

0.64 

2.00 

0.90 

SSZ (Sommer et ai, 1997) 

0 

1.20 

0.72 

2.20 

0.80 


In view of the temperature coupling between ( 1 ) andJ2), their closure requires know ledge of the varia- 
tion of tt 2 . This can be obtained by solving the modeled 9 2 and transport equations. Nine different sets of 
modeled equations are summarized in So and Speziale (1999) and they can be written in a general form as 



where P H = —»,(((<)(-) /9.v,). The different model constants adopted for pure and wall shear flows in the nine 
sets of modeled ^-equation are summarized in Table I, where abbreviations (shown in Table I) are used 
to denote the different models. These models are not completely independent of each other. For example, 
NLL and EL belong to one group. GJK. JM and SSZ another, AY, NK and YNT a third and CL is the fourth 
group. Members within a group only differ in the values of the constants chosen. In view of this, only one 
representative member of each group needs to be investigated. 


Numerical Solution 

Different combinations of the velocity governing equations given above will give rise to different Reynolds- 
stress closure models. For the sake of clarity, the following abbreviation are adopted to designate the various 
Reynolds-stress closures used for the calculations of turbulent buoyant shear flows. The 2-Eq designation 
is used to denote a two-equation model, which solves the trace of (I) and (10), the IP designation means 
adopting (3), (5) and (10), whije SSG means the use of (4), (5) and (10) in the equation set. It should be 
noted that, in the 2-Eq model, P is approximated by C,, (k 2 /e)S ii S lj with C fl = 0.09 and the constants asso- 
ciated with ( 10) are given by C ft = 1 .50, C l2 = 1 .90 and C f 3 — 1 .50. Other designations used are SSG/HGJ 
and SSG/SG. They signify the solution of the equation set (4), (6), (7) and (10) and the set (4), (6), (8), (9) 
and (11). respectively. Ot course, in all these closures, ( 1 ), (2) and ( 12)— ( 14) are also solved simultaneously. 
These closures differ in the modeling of the pressure-strain and dissipation rate tensors, and the pressure- 
scrambling vector. Comparisons of these models can. therefore, shed light on their ability or inability to 
predict equilibrium states of turbulent buoyant shear flows. 

The equations constitute a set of coupled initial value problems. They are solved using a Bulirsch-Stoer 
method (see Press et ai. 1986) with variable step size in order to achieve the desired accuracy with as few 
steps as possible. Calculations were carried out with the dimensionless equations. The physical equations 
were normalized assuming the length scale to be given by L, the time scale by I /S, the velocity scale by 
AC = SL and the temperature scale by AO = SL. where S = AU/dz. Therefore, the dimensionless time is 
given by r = St. It should be pointed out that even though the governing equations are not dependent on the 
molecular Prandtl number, Pr, the initial conditions are. The Pr dependence comes in through the specifi- 
cation of the initial which according to its definition involves the molecular thermal diffusivity. In view 
of this. Pr only affects the turbulence calculations in the initial stage of development but not in the long 
time period. The initial conditions are different for different flow' cases. Since the present study examines the 
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Table 2. The calculated values oiiP/e)^ and for homogeneous shear flows. 

Reynolds-stress closures 

2-Eq 



-j 

closure 

IP 

SSG 

SSG/HGJ 

SSG/SG 

(PM*, 

( Sk/F )oc 

1.886 

4.578 

1 .800 
4.927 

1 .886 
5.760 

1.886 

5.938 

1.901 

5.983 


approach to equilibrium turbulence of the cases considered by So et al. (1999), the same initial conditions 
specified there are used here. For most cases, unless specifically stated, Pr = 5 is assumed. 


Effect of Reynolds Stress Closures on the Equilibrium States of Pure Shear Flows 

The performance of different turbulence models in the predictions of homogeneous shear flows, homoge- 
neous buoyant flows and homogeneous buoyant shear flows has been critically examined (So et al., 1999) 
and the calculations were evaluated against DNS and experimental data. However, the approach to equilib- 
rium turbulence and the ability of these closure models to predict equilibrium states have not been examined. 
This paper proposes to examine the equilibrium states of stably stratified shear flows. Therefore, before pro- 
ceeding to analyze buoyant shear flows, the effect of the different Reynolds stress models on the prediction 
of equilibrium states for pure shear flows has to be examined first. 

For homogeneous shear flows, the approach to equilibrium turbulence is signified by dbij/dt = 0. This 
leads to P/e and Sk/e becoming constant as time goes to infinity (Jongen and Gatski, 1998). It may take 
a long time to reach this state, but large eddy simulation results (Rogers et al., 1986) show that such a state is 
achieved at equilibrium. Based on this condition and the governing k and e equations, the asymptotic values 
for P je and Sk/e can be determined. According to Jongen and Gatski (1998), the condition d(Sk/e)jdt = 0 
is valid at equilibrium. Then it follows that 

d (Sk/e) _ Sdk Skde_ { , 5) 

d/ e d / e 2 d t 

Substituting the trace of (I ) and (10) into (15) yields the asymptote. (P/e)oo = (G: - l)/(C f i - 1). An ex- 
plicit expression for Sk/e is possible under the assumption of P — C /t (k~ / e)S~ . This leads to (Sk/e)rx, 

l)/C M (C t ;i — l)] l/2 . Using the values of the constants specified for (10), (P/e) = 1886 and 
( Sk/e). x . =4.578 is deduced. This shows that production does not balance dissipation as equilibrium is 
approached. Instead, at equilibrium, P settles to a higher value than e. 

The next task is to evaluate the asymptotic values determined from the different Reynolds-stress closures. 
Only the model calculations from SSG and IP are considered, but they include the use of three different e,j 
models, i.e. those given by (5), (6) and (7), and (6), (8) and (9), respectively. The calculations are compared 
with the 2-Eq model results in Table 2. For the present calculations, the values of C F \ and Coadopted for the 
2-Eq. IP and SSG closures do not differ by much. Consequently, the asymptotic values of P/e determined 
from the model calculations cluster around 1.80. The correctness of these asymptotes for homogeneous 
shear flows remains to be verified. However, it should be pointed out that the experimental measurements 
of Tavoularis and Karnik (1989) indicate that (Sk/e) «, should fall between 5.0 and 6.0 and Speziale and 
Gatski (1997) also obtained a value for (Sk/e) close to 6.0. Based on this data, it appears that both the IP 
and the SSG results are fairly reasonable, irrespective of the dissipation rate tensor model assumed. All clo- 
sure models considered predict a (P/e) ^ value close to 1 .9 instead of 1 as a production-balance-dissipation 
model would imply. 


Equilibrium States of Buoyant Shear Flows 

For buoyant shear flows, the approach to equilibrium is further complicated by the presence of the thermal 
field. What additional parameters besides /ty, P/e and Sk/e should be used to characterize the equilibrium 
states? The DNS data of Gerz et al. ( 1 989) and the experimental measurements of Sirivat and Warhaft ( 1 983) 
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(a) T 



Figure 1, Variation of the components of h tJ with time assuming the JM model. 


could be interpreted to imply that the time rate of change of bjj and <y, should be zero as equilibrium is ap- 
proached. Furthermore, as the study of Zhao etal. (2001 ) showed, equilibrium buoyant turbulence as defined 
by dbij/dt = 0 and dcjj/dt = 0 will give rise to the state parameters P/e , Sk/e, Po/e e . R and G/e becom- 
ing constant. Therefore, this section investigates the approach to equilibrium turbulence as predicted by the 
different closure models. In anticipation of results to be presented later, only the velocity model calcula- 
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tions of IP, SSG, SSG/HGJ and SSG/SG, assuming an ^-equation model given either by JM or NK for 
a given Richardson number, Ri = f)g(d&/dz)/(dU/dz) 2 = 0.5, and Pr = 5. are shown in Figures 1-8. This 
calculated case is identical to the DNS data of Gerz cl al. (1989). The conclusions drawn from NK are also 
applicable to NLL and CL. Therefore, results obtained using NLL and CL and other members of the same 
group are not shown for the sake of brevity. 

The time variation of the components of b lt and q, calculated using a number of velocity models and 
the JM model for the % -equation are plotted in Figures 1 and 2. respectively. Only the results for the 
SSG, SSG/HGJ and SSG/SG are shown. It is obvious that the Reynolds stresses (Figure 1 ) have reached 
their respective constants as early as r = 100. As for the Reynolds heat fluxes, the stream component q j 
(Figure 2(a)) behaves like the Reynolds stresses, but the normal component qy (Figure 2(b)) reaches a con- 
stant value at r % 160. Thus, equilibrium is achieved only after r = 160. Once dbqjdt — 0 and dr/, /dr = 0 
have been achieved, R also reaches a constant value (Figure 3). The constant value reached is different for 
different velocity models. This is not unexpected, because the equilibrium values of (P/f)oc and ( Sk/e )r X> 
are different for different velocity models for pure shear flows as indicated in Table 1 . However, constant R 
does not necessarily imply that d(Sk/ s) /dt will go to zero also. It follows from the definition of R that when 
R and Sk/s are becoming constant, the following relation can be deduced: 

d(*/g) _ d(P/e„) Q (16) 

dr dr 



100 120 140 160 180 200 

(b) X 

Figure 2. Variation of the components of £/, with time assuming the JM model. 
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Therefore, plots of Sk/e and S0 2 /e versus x are also shown in Figures 4 and 5, respectively, to illustrate their 
behavior at large x. The results show that only the SSG/SG model yields constant Sk/e (Figure 4(a)) and 
sW/8 (Figure 5(a)) for r > 100. All other velocity models tested give an increasing Sk/e and S0 2 /e as r 
increases. Consequently, these results show' that the combination of models SSG/SG and JM is capable of 
predicting equilibrium turbulence for r > 160. Other velocity and thermal field models are not. 

If equilibrium is taken to be given by b (h q h /?, Sk/e and S0 2 /e all approaching constant, then the values 
of P/s and Py/ey can be derived as follows. Substituting the trace of (1 ), (10), (13) and (44) into (16), the 
following asymptotic relations are deduced for the equilibrium state parameters P / e and Py/ey. They are 

P C f 2 - 1 Qs - 1 / G \ (17) 

7 “ C F \ - 1 C F \ - 1 \e )' 

Py _ 2 C/3 R(C f 2 - 1 ) - 2 C f /5 R(C e 1 - 1 ) - (C d 4 - 2) (C ,1 - 1 ) _ 2C<hR / C^-\ \ /C \ 

77 ” (2 -C/I — 2C { pR)(C f \ - 1) 2-C/i -2CdiR \C F \ - 1 / V s ) 

( 18 ) 


so 2 

£g 



(a) 


T 



(b) t 

Figure 5. Variation of SO 2 ff-o with time assuming the JM model. 
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The expressions ( 1 7) and ( 1 8) relate P#/ e# and P /e to G/e and R and suggest that the equilibrium parame- 
ters Po/Sf) and P/e are constant if G/e is constant. It will be confirmed later that G/e is indeed constant at 
equilibrium, just as deduced by Zhao et al. ( 2CK) 1 ). 

In a buoyant flow, the normal heat flux would continue to oscillate until the collapse of turbulence (Still- 
inger et al.. 1983; Itsweire et al., 1986). Therefore, in principle, an equilibrium state could be reached only 
under the condition that G/e and other relevant parameters would become constant. This calculation shows 
that the SSG/SG and JM combination is able to predict P/e. Sk/e. Po/eh. R and G/e are constant at equi- 
librium, just as predicted by the analysis of Zhao et al. (2001). Unlike the case of pure shear flows, (17) and 
( 1 8) also depend on the constants specified in the modeled e- and ^-equation. In view of this, the constant 
values reached at equilibrium will depend on how correctly the e- and ^-equation can be modeled. 

Other combinations of velocity model and ^-equation model will not lead to a predicted equilibrium 
state. An example of one such combination (SSG and NK) is shown in Figures 6-8, where the time varia- 
tions of the components of bq and q,. and R are plotted, respectively. It can be seen that bq, q, and R fail to 
reach their asymptotic states even for very large r (however, only the period 1 00 < r < 200 is shown in these 
figures to illustrate the behavior). 


Effect of e„- Equation Modeling 

The above analysis on buoyant shear flows shows that at equilibrium bq. q : . P/e. Sk/e. P»/e H . R and G/e 
arc constant. Therefore, the effect of the ^/-equation on the prediction of the equilibrium state parameters 
is examined next. The calculations were mostly carried out using the 2-Eq and the Reynolds-stress models. 
In these calculations, the heat flux model invoked was that given by (12) with Cm set to zero. All e»- 
equation models listed in Table I w'ere investigated using the 2-Eq model. Only a few were examined using 
the Reynolds-stress models. The results for the 2-Eq model are shown in Figures 9-1 1, while those for the 
Reynolds-stressjnodels are plotted in Figures 3-5 and 8. From this point on. only the equilibrium values of 
R. Sk/e and S9~/e$ versus Ri are plotted. Their approach to equilibrium is not displayed because most have 
been shown in Figures 1-8. 

The 2-Eq model calculations were carried out specifying Pr = 5. Results show that not all modeled 
-equation can predict equilibrium turbulence. Of the four different groups of ^-equations investigated, 
only three predict equilibrium turbulence for ()</?/< 1. The value of R is plotted in Figure 9 for the six 
£ f y-equutionjnodels that yield equilibrium turbulence. Figures 10 and 1 1 show the equilibrium values for 
Sk/e and S0 2 /e„, respectively, over the same range of Ri. The R. Sk/e and S0 2 /e H thus calculated are differ- 
ent for different ^-equation model. In general, the values predicted by the CL model are greater than those 
given by the other models. Four models, GJK, JM, NK and YNT, yield very similar results, however. SSZ 
gives equilibrium values that are consistently higher than those obtained from GJK and JM. The correctness 
of these predictions will have to be verified by DNS data or experimental measurements. 

Not all Reynolds-stress closures can predict equilibrium turbulence with the four groups of ^-equations 
listed in Table 1 under the conditions dbq/<lt — 0 and dq,/dt = 0. Those ^-equations that cannot give equi- 
librium turbulence with the 2-Eq model also fail with the Reynolds-stress closures. Among the group of 
p^-equation models that can predict the approach to equilibrium turbulence, only the results of JM are shown 
because equilibrium is reached after a relatively short period of time. It took much longer time for GJK and 
SSZ to yield equilibrium turbulence. The results shown in Figures 3-5 are for the case where Pr — 5 and 
Ri = 0.5 (Gerz et al. 1989). All Reynolds stress models tested clearly show that R is constant (Figure 3). 
However, only SSG/SG give constant values for Sk/e and S0 2 jea (Figures 4(a) and 5(a)). The parameters 
Sk/e and S0 2 /e,, continue to vary with time long after R becomes constant for other Reynolds-stress models 
(Figures 4(b) and 5(b)). In other words, equilibrium has not been reached. It will be shown later that when 
R. Sk/e and S9-/eq become constant, so do P/e. P<i/eo and G/e. Other ^-equation models fail to give 
equilibrium prediction. A representative plot of this behavior is shown in Figure 8. 

According to the heat transfer studies made by the researchers of the nine different modeled in- 
equations. setting Cj2 = 0 gives better results for wall-bounded flows, while setting C,/| = 0 yields good 
agreement with the turbulence statistics of free shear flows. The other constants are found to have lit- 
tle effect on the calculations as long as they do not differ by much. Therefore, it appears that there is 
no one single set of constants that could give reasonably good results for both wall-bounded and free 
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Figure 6. Variation of the components of h tl with time for the SSG and NK models. 


shear flows. This is especially true when buoyancy-induced gravity waves (Sommer et al., 1997) and 
buoyant jets and plumes (Craft, 1991) are considered. In fact. Craft (1991) found that the constants 
c d 2 , C ( i4 and C d f, should be modified to make them dependent on R and the stress invariants. With 
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Figure 1 . Variation of the components of </,■ with time for the SSG and NK models. 



Figure 8. Variation of R with time for the SSG and NK models. 


this understanding, the inability of NLL, EL, AY, NK and YNT to predict equilibrium turbulence com- 
pared with JM is not surprising. As for CL, it assumes Cjy = 0. Therefore, it is reasonable to expect 
CL to perform differently from JM. The only difference between GJK, SSZ and JM is in the values 
assigned to the constants Cj2 to C,/.s. For example, GJK assumes a value of 1.80 for C,/ 5 , while JM 
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Figure 10. Variation of Skjt: with Ri at equilibrium for different ^-equation models. 
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Figure 11. Variation of S0 2 /ea with Ri at equilibrium for different ^-equation models. 

takes on 0.52. On the other hand, SSZ and JM invoke different constants for C d2 and C (/5 . Since all 
three models, GJK, JM and SSZ, are formulated for free turbulence, the present study shows that the 
choice of constants is important. Furthermore, the constants play a crucial role in the prediction of 
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figure 12. Variation ot R with Cm at equilibrium for different Ri. 
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the approach to equilibrium turbulence and of the asymptotic values of the state parameters in buoyant 
shear flows. 

Unlike the modeled ^-equation in the pure shear flow case, the modeled ^-equation affects the pre- 
dicted values of the state parameters and directly influences the prediction of the approach to equilibriim 
turbulence. 


Effect of Heat Flux Modeling 

Having examined the effect of the ^-equation model, the next task is to investigate the effect of <t> w model- 
ing. The above analysis is carried out with Cm = 0 in (12). In this section the effect of including this teim 
in the <t> w model on the prediction of equilibrium turbulence is studied. It should be pointed out that, in 
the calculations presented above, where Cm = 0 was specified m (12), both the IP and SSG closures using 
JM yield equilibrium turbulence where by, </,. R , Sk/e and S6 2 /e 0 are constant. Equilibrium turbulen:e 
was predicted for the range of Ri investigated. Therefore, only these two closures are investigated in detail 
in this section. The modeled equations solved are the same as before, except that the Cm term is retained 
in (12). The case studied is the DNS turbulent buoyant shear flow of Gerz et al. (1989) where Pr = 5.0 and 
Ri = 0.5. When the Cm term is added to (12), both IP and SSG closures predict equilibrium turbulence 
provided the value of Cm was modified to be 0.8. As for the value of Cm- it was found that a range af 
values would lead to the prediction of equilibrium turbulence. The plot of R versus Cm for the SSG clo- 
sure using JM is shown in Figure 12. The figure only shows the values of R for Ri =0.1, 0.5 and 1.0. This 
range of Cm values is not consistent with those recommended by Launder (1978). Therefore, these calct- 
lations tend to show that, when the pressure-scrambling term is modeled with the Cm term included, the 
IP or the SSG closure still predicts an approach to equilibrium turbulence for all values of Ri tested for 
a given set ot and Cm- However, the value of R calculated differs depending on the choice of con- 
stants. In other words, the addition of the Cm term to (12) has an appreciable effect on the prediction of 
equilibrium values. This effect is not desirable because it varies with the choice of Cm- Furthermore, it is 
not known which Cm would be more suitable for the developing stage of buoyant shear flows and for wall 
shear flows. 


Effect of Pressure-Strain Models 

The velocity closures are examined next and it is mainly concentrated on the effect of the pressure-strain 
models. Altogether two different pressure-strain models have been examined. These are the linear IP and 
the nonlinear SSG model. Calculations were carried out using the JM model for the ^-equation. The cal- 
culated equilibrium state parameters for different Ri are shown in Figures 13-16. All models (2-Eq. R 
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Figure 14. Variation of (a) P/e and (b) P w /fi« with Ri at equilibrium for different models. 


SSG. SSG/HGJ and SSG/SG) tested can predict constant bjj , </, and R. However, only the 2-Eq and 
SSG/SG closure are capable of predicting equilibrium behavior with Sk/s and S0 2 /eo also becoming con- 
stant (Figure 16). Together, these results show that IP, SSG and SSG/HGJ fail to predict the approach to 
equilibrium turbulence. In the case of IP and SSG, anisotropic modeling of the dissipation rate tensor is not 
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invoked. On the other hand, anisotropic modeling of the dissipation rate tensor is assumed in SSG/HGJ and 
SSG/SG. The differences are in the models assumed in djj and in the modeled e-equation. In SSG/HGJ the 
anisotropy is expressed as a function of bjj, while SSG/SG represents djj in terms of Sjj and Wjj and solved 
a different f -equation. The present results seem to imply that either the form used in SSG/SG to account for 
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dissipation anisotropy or the ^-equation assumed is more appropriate. It will be shown in the next section 
that a simple modification of the e-equation is sufficient to render IP and SSG valid in the prediction of the 
approach to equilibrium turbulence. This may not be surprising since the equilibrium analysis did not need 
to account for dissipation rate anisotropy. 

Asymptotic values of P/e and Ph/sh thus determined for different Hi are shown in Figure 1 4. In the case 
of p/e. the 2-Eq prediction is of opposite trend to that given by the SSG-type closures. Instead of a decreas- 
ing trend for increasing Ri, the 2-Eq result shows an increasing behavior with Ri. According to the SSG/SG 
calculations, the heat flux contribution is quite small and is positive, while the shear component is large and 
is negative. Therefore, the trend displays by the 2-Eq closure could not be correct. As for the behavior of 
Ph/£h (Figure 14(b)). again the SSG-type closures yield a decreasing trend similar to that of P/e. However, 
the value reached at Ri = 1 is much smaller than that given by the 2-Eq closure. This implies that stratifica- 
tion suppresses the production of 6- more strongly than its dissipation. On the other hand, the 2-Eq closuie 
predicts that the production of 6 2 and its dissipation are roughly equal for Ri > 0.5. The behavior ot G/e 
with Ri is shown in Figure 1 5. Just as before, the values predicted by the SSG-type closures are quite similar, 
while those given by the 2-Eq closure essentially decrease linearly with Ri. 


Effect of Anisotropic Dissipation Rate Modeling 


The next task is to examine the effect of anisotropic dissipation modeling on the calculation ot equilibrium 
states. Three different dissipation models are considered. The first is the isotropic model (5). the second is 
that given by (6) and (7), while the third is provided by (6), (8) and (9). The IP and SSG closures are used 
to carry out the calculation. In all these calculations, the JM model is invoked tor the ^-equation and ( I -) is 
solved with Cyn = 0. Again, the case calculated has the same initial conditions as those given by Gerz et al. 
(1989). 

From the above analysis, it is known that different dissipation rate models give rise to different predic- 
tions of the equilibrium values. The reason why the SSG/SG closure can predict equilibrium turbulence 
while the IP. SSG and SSG/HGJ closures cannot is, perhaps, due to a difference in the modeled e-equation 
invoked by the closure models. The SSG/SG closure assumes (11) while IP, SSG and SSG/HGJ invoke ( 10). 
The difference between these two equations is the additional d,j .S',, term in (1 1). In conjunction with (8) and 
(9), the leading term of iljjSjj is given by a constant times kSjjSjj. This term is similar to a vortex stretch- 
ing term introduced by Bernard and Spe/.iale (1992) to drive the flow to a production-balance-dissipation 
equilibrium with bounded energy states. Thus (11) could be interpreted as a first attempt to account for 
anisotropic dissipation, yet also serves to strain the vorticity by mean shear. The question then is could such 
a term in the modeled e-equation be sufficient to remedy the shortcoming of the IP and SSG closures without 
having to invoke an anisotropic dissipation rate model in ( 1 ) such as given by (8) and (9)? An attempt is made 
to answer this question here. 

An e-equation with an additional term to account for vortex stretching has been proposed by Bernard and 
Speziale (1992). It is derived from the equation governing the behavior of enstrophy or — co/co/ in isotropic 
turbulence and the identity e = vor. Thus derived, the additional term depends on the fluid kinematic viscos- 
ity, v. and the modified e-equation for turbulent buoyant shear flows becomes 


de 

dr 


e - 7 

= C ftT P+— r - r 
k 3/15 xA 1 


c* 


-C f -2y . 


(19) 


The model constants take on the following values: C t \ - 1 .44, C f 2 = I -83, C* 3 - 0.01 and C t: y — 1 .44. The 
vortex-stretching term renders the solution a function of the initial turbulent Reynolds nuinbei Re (So et ul., 
1999). This is obviously not the case, even for pure shear flows. An alternative to this proposal is given 
by (11). The additional term in (11) could be simplified to give a term quite similar to the one proposed 
by Bernard and Speziale ( 1992) and yet not dependent on the initial turbulent Re. Since the leading term 
in edjjSjj is a constant times kSySy and k 2 SijSij/e 2 is a constant for homogeneous shear flows, the term 
k 2 (SijSjj) 2 could be related to e 2 SjjSjj. Thus modified, the e-equation (11) could be written as 

2 

— = C y] jP + C,4 s/2( Si, Sjj ) 1 /2 e - C, 2 Y - C F y j . 
dr k k k 


(20) 
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Figure 17. Variation of R with Ri at equilibrium for different r-cquations. 
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where C F 4 — 0.042 is a model constant that yields the best predictions for pure shear flows. 

In the following, AD is used to denote the c-equation (20) and the SSG closure that solves (20) is 
designated as SSG/AD. In this closure the dissipation rate tensor is still modeled by (5). Thus, SSG/ AD 
represents a first attempt to approximate the effect of anisotropic dissipation. If the results come out to be 
similar to those given by SSG/SG, then there is no need to solve (6). (8) and (9). The performance of (20) 



Turbulence Modeling Effects on the Prediction of Equilibrium States of Buoyant Shear Flows 


419 


G 

£ 


0 

- 0.1 
- 0.2 
- 0.3 
- 0.4 
- 0.5 
- 0.6 

0 0.2 0.4 0.6 0.8 1 


JM model 

• SSG 

A SSG/AD 
▼ SSG/BS 
■ IP/AD 

♦ SSG/SG 


Ri 

Figure 19. Variation of G/f with Ri at equilibrium for different /--equations. 
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Figure 20. Variation of (a) Sk/e and (b) S# 2 /e ( , with Ri at equilibrium for different f -equations. 

has been validated against the measured decay ot k/ko with time (Tavoularis and Karnik, 1989) by So et ai 
(1999), where fc 0 is the initial value of k. Other closures used to calculate the same flow include 2-Eq, IP, 
SSG, SSG/HGJ and SSG/BS. The SSG/BS closure is the same as SSG but solves the ^-equation (19) in- 
stead of (10). The SSG/AD and SSG/BS closures were found to yield good agreement with experimental 
measurements. All other closure models did not. However, the SSG/BS closure sutlers from the shortcoming 
of being dependent on the initial Re. This cannot be the case for long time. 
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With a new proposal to account for anisotropic dissipation, the ability of the SSG model to predict equi- 
librium states of turbulent buoyant shear flows is greatly improved. If (20) instead of (10) is used with the IP 
and SSG models, the calculations show that they can predict the approach to equilibrium turbulence much 
like the SSG/SG closure, i.e. at equilibrium P/e, P H /e {i , R, G/e, Sk/e and SV 2 /e tt all approach constant. 
The results of these calculations are shown in Figures 1 7-20. Altogether, five closure calculations are show n, 
they are SSG. SSG/AD, SSG/SG, IP/AD and SSG/BS. The variations of R, P/e, P«/e 0 , G/e , Sk/e and 
SO- j eh with Ri are plotted in Figures 17, 18(a), (b), 19, and 20(a), (b), respectively. From these results, it 
is obvious that SSG/BS cannot predict equilibrium turbulence with R , Sk/e and S9 2 /s ( ) becoming constant 
simultaneously, even though P /e, Pv/eh and G/e also achieve constant behavior. It appears that the R pre- 
diction ot SSG/AD is very close to_that of SSG/SG. This is also true for 1P/AD. On the other hand, the 
SSG/AD predictions of Sk/e and S9 2 /e H are higher than those given by SSG/SG by a factor of 2. The pre- 
dicted behavior and magnitudes of P/e, P^/e 0 and G/e given by SSG/SG, SSG/AD and IP/AD are very 
similar. It can be seen that all closures tested give a decreasing trend for P/e and Po/sn with Ri , but an 
increasing trend with Ri for R. The predictions, therefore, imply that the production of A and 9 2 decreases 
laster than their dissipation as Ri increases. These results tend to suggest that the modeling of e t j has an 
effect on the prediction of equilibrium turbulence. It is important to account for anisotropic dissipation, how- 
ever, it is not necessary to model the effect in £,■/ directly; it could be approximately accounted for in the 
e -equation. 


Conclusions 

This investigation ot the effects ot turbulence models on the prediction of equilibrium states of buoyant shear 

flow' leads to the following conclusions: 

(a) The equilibrium turbulence state defined by dhjj/dt = 0 and dqifdt = 0 leads to R, P/e, fy/eu, G/e 
and Sk/e becoming constant at equilibrium, consistent with the analysis of Zhao et al. (2001 ). 

(b) Among the different modeled ^-equations considered, only three groups could lead to the prediction 
of equilibrium turbulence when a 2-Eq closure is used. If the velocity held is closed by SSG/SG. only 
one group of models could give rise to equilibrium turbulence. In this group only the JM model would 
predict equilibrium turbulence after a relatively short period of time. 

(c) All other closures, IP, SSG, SSG/SG and SSG/HGJ, even in conjunction with JM, could not predict the 
approach to equilibrium turbulence. 

(d) It is found that variations in the modeling of the pressure-strain tensor and the pressure-scrambling 
vector have little or no effect on the prediction of the approach to equilibrium turbulence. The actual 
predictions of the equilibrium values may differ slightly though. 

(e) The SSG/AD and IP/AD closures can predict the approach to equilibrium turbulence just like that of 
SSG/SG. This suggests that a simple accounting for anisotropic dissipation behavior could be accom- 
plished through a suitable modification of the ^-equation. 
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